Self-Organized Critical Random Boolean Networks 
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Standard Random Boolean Networks display an order-disorder phase transition. We add to 
the standard Random Boolean Networks a disconnection rule which couples the control and order 
parameters. By this way, the system is driven to the critical line transition. Under the influence of 
perturbations the system point out self-organized critical behavior. Several numerical simulations 
have been done and compared with a proposed analytical treatment. 
PACS numbers: 05.65.+b, 87.23. Ge, 87.23.Kg 



I. INTRODUCTION 

Random Boolean Networks (RBN) were proposed Jl]] 
as discrete genetic networks models. The network is com- 
posed of nodes (automata). The state of an automaton 
represents the state of a gene with two possible values 
(on-off). This state is the output of a Boolean function, 
which has as inputs the output activity of some other 
genes. The connectivity of the system and the bias used 
for the Boolean functions are relevant parameters in or- 
der to statistically determine the network dynamics. If 
the system has a high connectivity and a low bias, the dy- 
namics of the automata is disordered; it seems that there 
is no correlation between the gene switching on and off. 
On the other hand, the dynamics is ordered if a low con- 
nectivity and a high bias are used. Only the ordered dy- 
namics have biological sense. S. Kauffman says that this 
parametric region offers "/ree orrfer" and it seems that 
natural selection works where the order exist previously 
H . In the critical region, which is the boundary between 
ordered and disordered phase, there are some scaling re- 
lations that have been the subject of recent works ||-|| ■ 
S. Kauffman points out that genetic networks evolve to 
the boundary between order and disorder. In this re- 
gion there is more diversity of patterns for activity and 
greater possibilities for complex evolution (anti-chaos hy- 
pothesis) ||. 

There is great interest in the evolution of topology in 
networks MM. Several schemes have been proposed for 
RBN evolution. For instance, in j^] a RBN evolves from 
an initial K = 1 mean connectivity (using a random ini- 
tial condition) to an attractor. At this point a copy of 
the network is made which has a connection randomly 
removed or /and added. The new network also reaches 
an attractor. If the same attractor has been obtained, 
the network is maintained. Otherwise, the previous one 
is restored. During the RBN evolution, it appears some 
stasi periods and punctuations like in real evolutionary 
processes. 

In the performance of an automaton is defined as 
the number of steps that it is in majority during a given 
set of steps. An automaton in majority means that it 
has the same state as the majority of the automata. The 
automaton with the highest performance is replaced by 
another automaton with a new random Boolean function 



and the process is repeated. The authors show how the 
genetic network is able to modify its bias in order to reach 
the critical region. 

In jllj is presented a method which is able to lead 
the connectivity of binary neural networks to their crit- 
ical value. These neural networks have a phase transi- 
tion similar to RBN jl2],[l3|. Using a fixed connectivity 
initialization the network evolves to an attractor. If an 
automaton does not change its state in the attractor, a 
connection is added to it. Otherwise, if it changes then 
a connection is removed. The different networks (with 
different initial connectivities) reach and remain close to 
the critical region. 

From a more general perspective, some works jl4|-|l8]l 
have shown the relation between critical and self- 
organized critical phenomena. In particular, D. Sornette 
et al. yl| point out a heuristic method that transforms a 
system with a critical transition (Ising or bond percola- 
tion) into a self-organized critical system. They propose 
to add some kind of mechanism (a feedback between the 
order and control parameters) which slowly drives the 
system towards the critical point. 

RBN are a classical example where complex global be- 
havior emerges from local simple rules. They exhibit a 
phase transition similar to the Ising model or bond per- 
colation. Forgetting its initial biological inspiration, our 
goal is to develop RBN which spontaneously evolve to- 
wards a global critical stationary state. For this purpose 
we use a disconnection rule which induces a feedback 
between the control and order parameters. These net- 
works reach a critical state without changing externally 
the control parameter. These systems show character- 
istics related to self-organized criticality. Our evolution 
method is distinct from earlier works in the following: 
we introduce a well-defined coupling between order and 
control parameters as in p| , but the method is able to 
stabilize (the individuals connectivities become constant 
in time) the networks along the critical curve. 

In section II we give an introduction to RBN and we 
present the disconnection local rule which gives to the 
RBN a SOC behavior. In section III statistical results 
of the evolution of the self-organized RBN are presented. 
In section IV we show its SOC behavior in response to 
external perturbations. In section V we apply an analyt- 
ical treatment to the model and discuss the results, and 
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finally, in section VI we make a summary and point out 
for future works. 



2. Minimal connectivity: if Ki(t) = 2 then 



II. HOW DOES ONE MAKE A CRITICAL 
SELF-ORGANIZED RBN 



Ki(t + l) = Ki(t). 



(7) 



A RBN is a discrete dynamical system composed of N 
automata. Each automaton is a Boolean variable with 
two possible states: {0, 1}, such that 



F: {0,1}^ ^{0,1} 



/v 



(1) 



where F = — ,/»,—, /iv) an d each fi is a Boolean 
function of Ki inputs (the automaton i is connected to Ki 
automata randomly chosen from the set of N automata) : 



fi : {0,1}*' .-{0,1}. 



(2) 



An automaton state x\ S {0, 1} is updated using its cor- 
responding Boolean function: 



r.t+1 



fi (S3 jj , X j 2 , . .. , %i Ki ) • 



(3) 



Each fi is represented as a look-table of Ki inputs. Ini- 
tially, Ki neighbors and a look-table of bias p are as- 
signed to each automaton. In order to generate such a 
look-table, the value 1 is assigned to an output with a 
probability p and with a probability 1 — p. When the 
neighborhood and the functions are established, they are 
maintained (quenched). 

We randomly initialize the states of the automata (ini- 
tial condition of the RBN). The automata are updated 
synchronously using its corresponding Boolean functions: 



F(x* 



(4) 



These RBN exhibit a second order phase transition 
JUSf . The control parameters (K and p) determine two 
regions: a frozen phase for K < l/2p(l — p) and a dis- 
ordered phase for K > l/2p(l — p). Thus, the critical 
boundary is represented as follows: 



K c {p) = 



1 



2p(l-p)' 



(5) 



The above description corresponds to a classical RBN. 
We incorporate a disconnection rule to the system lead- 
ing the RBN to a stationary critical state with connec- 
tivity K c (p) (see equation ||). The rule is also applied 
synchronously to each automaton i which has a local con- 
nectivity, Ki(t), in such a way that: 

1. Disconnection threshold: if Ki(t) > 2 then 

K,(t) 

Ki(t + 1) = Ki(t) - 1, if x i, < K i(t) - 1 

j=i 

Ki(t + 1) = Ki(t), otherwise. (6) 



The rule is inspired in the Bak-Tang-Wiesenfeld 
(BTW) model pfj. As in the sandpile-like model we 
have incorporated a threshold. If the number of 0's in 
the input of a Boolean function is greater than one and 
the local connectivity of its corresponding automaton is 
Ki{t) > 2, a connection of this automaton (randomly 
chosen) will be cut. Then, a new Boolean function with 
connectivity K%{£) — 1 will be assigned. 

The disconnection mechanism is illustrated in Figure 
[j]. The automaton i has a connectivity Ki(t) = 3 and an 
input vector composed by two 0's and one 1. Hence, a 
random connection is cut (Ki(t + 1) = 2) and another 
Boolean function with bias p = 0.5 is assigned. A new 
state is then computed. 

In a sandpile-like system there is a decrease in the 
mass, energy, tension, etc. at the boundaries p5] , pl| . 
In the proposed RBN, there is a decrease in the number 
of connections. Our system is a random graph with no 
boundaries. But, there is a diffusion effect which is able 
to cause an avalanche of disconnections. If one automa- 
ton connection is cut (at t), its state can change from 1 to 
(at t + 1). The last change may cause a disconnection 
(at t + 1) in any of the automata that are connected to 
the first one, and so on until a stationary state is reached. 
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FIG. 1. Disconnection of an automaton i with Ki(t) — 3. 
The automaton receives more than one 0. At t + 1/2 one of 
the automaton connections is randomly cut and its Boolean 
function is changed (by a look-table with bias p = 0.5). The 
new Boolean function obtains a new state for the automaton 
in t + 1. 
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FIG. 2. Average evolution of the mean connectivity K(t), 
where t are the iterations. We have considered 1000 RBN 
(N = 10000): 100 different RBN with 10 initial conditions 
(50% of the automata with state 1 and 50% with state in 
t = 0) for each value of the bias p. All the RBN have a initial 
K(t = 0) = 10 mean connectivity. 
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FIG. 3. (a) The continuous line represents the critical 
curve obtained from the theory (equation ^). It is a bound- 
ary between the ordered and disordered phases of the RBN. 
(b) The fill down-triangles represent three instances of the 
RBN evolution for p = 0.50, p = 0.70, and p = 0.85. Each 
triangle represents the mean connectivity K (t) in consecutive 
steps, (c) The circles represent the stationary state at the 
end of the RBN evolution (K(t = 0) = 10) using different 
values for the bias p (ranging from p = 0.5 to p = 0.9 in p in- 
crements of 0.01). Each circle represents the average of 1000 
RBN (N = 10000) during 1000 time-steps. 



III. RESULTS 



The average of the evolution of the mean connectivity 
K(t) of 1000 RBN, each one with N = 10000 automata, 
is represented in Figure Each RBN diminishes its con- 
nectivity due to equations || and |t]. A fixed bias p has 
been used for every RBN: 100 different RBN and for each 
one 10 initial conditions of the automata states (50% of 
l's and 50% of 0's approximately). All the RBN au- 
tomata have a connectivity Ki(t = 0) = 10. The bias 
p ranges from 0.5 to 0.9 in steps of 0.01. The station- 
ary state (when no more connections are cut) is quickly 
reached in about 50 iterations. The mean connectivity 
K(t) stabilizes close to the transition curve between or- 
der and disorder (K c (t)). This curve is determined by 
equation |^ and it is represented as a continuous line in 
Figure H In the same figure, the circles are the stationary 
result of the averaged mean connectivity (corresponding 
to section p versus K at t = 100, in Figure ||). It can be 
observed that the RBN stabilizes at the transition curve. 
In the same figure, three instances of the evolution of 
three RBN (p = 0.50, p = 0.70 and p = 0.85) have been 
represented. The fill-down triangles represent the values 
of the RBN mean connectivity evolving towards equilib- 
rium. In addition the patterns evolution of the above 
three examples are represented in Figure 0. 



The self-overlap a(t) is the unitary percentage of au- 
tomata with the same value in t — 1 and in t. The station- 
ary self-overlap, a(t — > oo) = a*, is an order parameter 
for the RBN transition p2|. Therefore in the disordered 
state we have a* < 1 and in the ordered state the self- 
overlap is given by a* = 1 . In Figure ^| we have repre- 
sented the self-overlap evolution (averages are calculated 
as in Figure ||). The initial self-overlap is a(t = 1) = 0.5 
due to the random state initialization. It can be observed 
that the self-overlap grows quickly reaching the station- 
ary value a* = 1 for high values of the bias. This is to be 
expected because in the critical boundary with K = K c , 
the self-overlap reaches the value a(t) = 1 for the first 
time. On the other hand, for low values of the bias, the 
self-overlap is stable at a value lower than 1. We will 
show in section V that this is an RBN size effect. 

In order to compare the different velocities to approach 
the stability in the self-overlap and the connectivity, the 
self-overlap of the connectivities q(t) has been repre- 
sented in Figure ^. The quantity q(t) is defined as the 
unitary percentage of automata which have the same con- 
nectivity at times t — 1 and t. As it can be observed, 
the connectivity is stabilized, approximately, in 50 time- 
steps. On the other hand, the self-overlap needs about 
100 time-steps to stabilize. 

In short: the evolution of an RBN with dynamics given 
by equations [||?] is such that the RBN is driven from the 
disordered phase to the critical boundary (equation ||) 
with mean connectivity K c (p) and self-overlap awl. 

In the next Section we study the nature of the RBN 
equilibrium state, and whether it is meta-stable or not. 
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In order to do so we analyse the response of the system 
to external perturbations. 




FIG. 4. These figures show the spatio-temporal automata 
activity for the evolution of networks built of N = 100 au- 
tomata. Time runs from left to right. The number of it- 
erations is represented in the figure by t. The value of the 
automaton state (vertical axis) is represented as black if it 
is one, and white otherwise. From top to down: RBN with 
p = 0.50, p = 0.70, and p = 0.85. 




0.5 



FIG. 5. Evolution of the average self-overlap a(t) for the 
simulations of Figure (t is the number of iterations). 



We now evolve the RBN using the disconnection rule 
(equations || and 0) . We perturb the RBN making use of 
sandpile-like methods [ po|j2l| . For this purpose we start 
with a relaxed RBN and a fixed bias p. In order to per- 
turb the RBN, we randomly choose an automaton i, and 
add to it a connection (i.e., its connectivity changes from 
Ki to Ki + 1). A new Boolean function fi with Ki + 1 
inputs is assigned to i using the bias p, so that one ob- 
tains a new state for the automaton i. If the automaton 
state is maintained, then there is only a small increase 
in the RBN mean connectivity. If the automaton state 
changes, it is possible to cause a disconnection, origi- 
nating an avalanche via branching. When the avalanche 
stops, the RBN is minimally perturbed again, and so on. 

In order to characterize the avalanches, we have mea- 
sured two different variables: the time T needed for the 
net to reach a stationary state and the total number of 
disconnections B during T. In Figure the histogram 
S(T) of avalanche times is represented after 5 x 10 6 
perturbations. For this purpose we have used a RBN 
with p = 0.65 and different sizes: N = 100 (triangles), 
N = 1000 (unfilled circles), N = 10000 (filled circles), 
and N = 100000 (diamonds). The dotted line is a visual 
guide to see the size effect. It is easy to see that when 
N — > oo, the histograms tend to a power law. A fit with 
the first points gives S(T) ~ T z with z = -2.9 ± 0.1. 
Figure || shows the histogram S(B) of disconnections us- 
ing the same numerical simulations as before. We obtain 
again a possible power law described by S(B) ~ B D with 
D = —2.2 ± 0.2. In both figures, insets in linear-log form 
are showed in order to exclude posible exponential fits. 
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FIG. 7. Log-log histogram (S(T)) for T (iterations needed 
for the system to reach a stationary state) during 5 x 10 6 per- 
turbations in a RBN. p = 0.65 for N = 100000 (diamonds), 
N = 10000 (filled circles), N = 1000 (unfilled circles), and 
N = 100 (triangles). Inset: the same figure but in linear- log 
form. 
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FIG. 8. Log-log histogram for the sizes of the avalanches 
B. B is the number of disconnections for each perturbation 
and S(B) is the histogram. For this purpose it has been used 
the same simulations of Figure (?]. Inset: the same figure but 
in linear-log form. 

In the inset of Figure ^ we show the temporal evolu- 
tion of the mean connectivity K(t) for a perturbed RBN 
(N = 10000, p = 0.65). In Figure ^ we represent the 
power spectrum S(f) (log-log) obtained by averaging 305 
temporal series of 16384 time-steps. A power law de- 
scribed by S(f) ~ / with = -1.92 ± 0.09 is in good 
agreement with the data. 

The numerical results do not seem definitive. The 
power law regime for T and B lasts almost for a decade, 
but there is a clear finite size effect. There are strong 
computational restrictions for working with bigger RBN 
because the size of the avalanches is limited by the num- 
ber of RBN automata. Nevertheless, in the S(f) his- 
togram there is scaling for more than two decades. In 



conclusion, we think that the results point out a SOC 
behavior. 
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FIG. 9. Power spectrum (log-log) of the mean connectivity 
evolution K(t) of a RBN (p = 0.65 and TV = 10000). It is the 
average of 305 temporal series of 16384 steps. The slope of 
the dotted line is <j> — — 1.92±0.09. Inset: a particular case of 
the evolution of the mean connectivity K (t) for the average 
series (t is the number of iterations). 



V. ANALYTIC TREATMENT 

At this point we ask ourselves: why are the evolution 
rules able to drive an RBN to the critical boundary be- 
tween the order and disorder phases? As the bias has 
been already fixed, there can only be variations in the 
mean connectivity K(t) (acting as order parameter) and 
the self-overlap a(t) (control parameter). In this section 
we show that there is a feed-back mechanism between 
both parameters that leads the system to a SOC behav- 
for @. 

As we mentioned before, the self-overlap a(t) is the 
unitary percent of automata which has the same state in 
t + 1 and t. If the mean connectivity K(t) is known, the 
value of the self-overlap at time t + 1 can be determined 
by the following equation: 

a(t+l)=a K V>(t)+V(l-a K ®(t)). (8) 

If we interpret a(t) as the probability for an arbitrary 
automaton to remain in the same state at both t—1 and 
t, the term a K ^\t) gives the probability that all the in- 
puts of a given automaton are the same for t— 1 to t. It is 
clear that 1 — a K ^'(t) is the probability that at least one 
of the inputs will be different at t and t—1. In that case, 
there is still a probability that this particular automaton 
remains in the same state at t and t + 1 by chance. This 
probability is given by V = p 2 + (1 - V? H- In a stan- 
dard RBN the mean connectivity K does not have an 
explicit dependence on time. If one includes the discon- 
nection rule described above, the mean connectivity of 
the system evolves with time. If Pfc(£) is the probability 
of an automaton to have connectivity k, then: 



5 



10 



(9) 



k=2 



where the maximum connectivity value has been taken 
as 10, without loss of generality. The connectivity distri- 
bution evolves according to the following system of equa- 
tions: 

P w (t + 1) = [$ 10 + (1 - $ 10 )a 10 (i)]Pio(i) (10a) 

Pk(t+1) = (1 - *fc+i)(l - a k+1 (t))P k+1 (t) 

+ [$ fe + (1 - 1> k )a k (t)}P k (t) (10b) 

P 2 (t + 1) = (1 - * 3 )(1 - a 3 (t))P 3 (t) + P 2 (t), (10c) 

where k = 9, 8, 3. The variable Q k — P k + kp k ~ 1 (l —p) 
is the automaton probability to maintain its connectivity 
by chance. Equation 10a describes the loss of connec- 



tions for all the automata wi th co nnectivity 10 at instant 
t. T he rest of the equations |l0b| , except for the last one 
10c , have two contributions: the first one represents the 
creation of automata with connectivity k that previously 
had connectivity k + 1, and the second one describes the 
automata that maintain their same value of connectivity 
k. The last equation 10c describes the growth of au- 
tomata population with connectivity 2. The evolution of 
the critical self-organized RBN is descr ibed by the cou- 
pled system of equations || and 10a - 10c . 

One can see than a* = 1 is the only possible value that 
makes the probability distribution become stationary. In 
equation (g), a* — 1 is a fixed point when K*(l — V) < 1, 
where K* = K(t — > oo) is the asymptotic connectivity, 
that is, whenever the stationary mean connectivity sat- 
isfies: 



K* < 



1 



(11) 



whose boundary is the critical c urve dep icted in Figure 
||. From the system of equations lOa- lOoj it can be easily 
deduced that: 



K(t+l) = K(t)-(l-P 2 (t)) 

10 

■X)[** + (i-**)«*(*)]^(<). 

fc=3 



(12) 



so that the mean connectivity of the system always de- 
creases during its dynamical evolution. Therefore, if the 
initial condition satisfies K(t — 0) > l/(2p(l —p)), it will 
slowly fall towards the critical condition. The first time 
a* = 1 is when the system reaches the critical curve and 
it, therefore, stabilizes. 
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FIG. 10. The theoretical evolution of the self-overlap 
in a RBN. Where t is the number of iterations, p — 0.5, 
a(t — 0) — 0.5, and Kit — 0) = 10. All the automata have 
connectivity 10, i.e., Pio(t = 0) = 1.0 (continuous line). The 
circles correspond to the numerical simulation for N = 10000, 
the triangles for iV = 1000 and the squares for iV = 100. In 
the inset, the filled circles represent the average stationary 
values of the self-overlap reached by the RBN with different 
sizes, and logarithmically spaced out. 

We can study more rigorously the evolution of the sys- 
tem by the classical analysis of linear stability. We per- 
turb the fixed point ft* = (a* = 1, Pf , P*, P 2 *), 

with K* = Sfc=2-Ffc^' and compute the Jacobi matrix 
of the system at ft*: 



L(ft*) 



K*{l-V) 
10(1 - $io)Pi* 







fe(l - 
-(fc + l)(l 

-3(1- 



-^+i)p; 



k+l 



4> 3 )P 3 * 







1 : : 

: 1 
1/ 

(13) 



We can calculate the characteristic polynomial by means 
of the determinantal condition: 



Thus: 



|L(ft*) - All = 0. (14) 



P(\) = (K*{l-V)-\)(l-\f = 0, (15) 



and then the associated eigen- values are Ai = K*(l — V) 
and the (9 times) degenerated eigen- value \ 2 = 1. This 
means that for Ai < 1 the critical point will indeed be lin- 
early stable. One can see that the characteristic polyno- 
mial is independent of the stationary distribution {P k }- 
Obviously this distribution depends on the initial condi- 
tions, but fulfills the condition K* ^* 1 
the condition |TT| . 
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FIG. 11. Theoretical evolution (continuous line) and nu- 
merical simulation (circles) of a(t) for RBN (N = 10000). We 
have performed 1000 averaged simulations with four different 
values for the bias p (t is the number of iterations). 



that the stationary state of the simulations approaches 
the theoretical one as N increases. In order to perform 
a scaling size effect study we have represented the sta- 
tionary states (1 — a*) for different sizes (see the inset 
in Figure [h]). The fitting shows that (1 - a*) ~ (1/N) a 
with a — 0.27 ± 0.05, in such a way that in the ther- 
modynamic limit (N —> oo) we get a* = 1, as is to be 
expected. 

In a similar way, the four graphs in Figure [ll] show the 
theoretical evolution (continuous line) and the numerical 
simulation (circles) of the overlap for RBN (N = 10000). 
For each graph we have used 1000 averaged simulations 
and four different values for the bias p. In Figure ^ we 
have plotted the evolution of the theoretical K[t) mean 
connectivity (lines) and simulations with RBN (symbols) 
for the same bias values of Figure 11. Finally, in Figure 
|l3| we show {P£} the stationary distributions as func- 
tions of the bias p: theoretical (top one ) and simulation 
(bottom). 
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FIG. 12. Evolution of the mean connectivity K(t) (t is the 
number of iterations): the lines correspond to the theoretical 
results and the symbols to the averaged simulations. We have 
considered the same cases as in Figure |llj 

In order to check this theoretical analysis we have per- 
formed several numerical simulations. In Figure [l^ the 
theoretical evolution of the self-overlap a(t) has been rep- 
resented with continuous line. We have used the following 
values: a(t = 0) = 0.5, p = 0.5, and P k (t = 0) = 0.0 for 
all k but Pio = 1.0, that is, K(t = 0) = 10. It can be 
seen that the self-overlap converges asymptotically to 1. 
With the same parameters and initial conditions we have 
performed our set of numerical simulations. In the same 
figure we have also plotted the average evolution of the 
self-overlap for 1000 RBN of different sizes: TV = 10000 
(circles), N — 1000 (triangles), and TV = 100 (squares). 
As can be seen from the plot, the first steps of the simula- 
tion agree with the theoretical curve. It can be observed 




FIG. 13. Stationary distribution of the connectivity P*(k) 
with k = 2,3, 10. In the bottom graph each distribution 
has been calculated averaging over 1000 RBN. In the top one 
we show the theoretical stationary distribution of connectivi- 
ties. 
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From the previous figures it can be clearly observed 
that the computer RBN simulations agree with our the- 
oretical analysis, in spite of the size effect, that is due to 
computational restrictions. 



VI. SUMMARY 

RBN are classic complex systems which show that 
global order is able to emerge from local rules. They 
exhibit a phase order-disorder transition modulated by 
the values of their connectivity and bias. Some authors 
have already proposed mechanisms for RBN evolution 
[ gjUj . We have, following the Sornette et al. criterium 
[ njTcoupled the control parameter (connectivity for a 
fixed bias) and the order parameter (self-overlap) in order 
to convert a critical system into a self-organized critical 
one. From the numerical simulations it can be deduced 
that the system evolves from the disordered phase to the 
critical transition curve independently of the fixed bias. 
Moreover, the response to external perturbations points 
to a SOC behavior. 

The theoretical analysis fits with the numerical simu- 
lations and it is able to predict RBN behaviors to govern 
variations. We tested other disconnection rules (similar 
to ^| and 0): for example, we change the breaking thresh- 
old of a connection. If there is a disconnection when no 
0- inputs arrive to the automata, the system stabilizes un- 
der the critical curve. If there is a disconnection when 
there are more than one 0-inputs in the automata, the 
system stabilizes over the critical curve. In the theoreti- 
cal analysis, this modifications imply that the probabili- 
ties &k change. In the first case the probabilities increase 
and in the second case they decrease. The values of these 
probabilities control the rate at which the mean connec- 
tivity decreases. In the first case, the rate is greater and 
the self-overlap reaches the value I deep inside the or- 
dered phase. In the second case, the rate is very small 
and the net stabilizes in the disordered phase with self- 
overlap lower than 1. There is a range for the values of 
the $fe functions that gives the correct rate. The chosen 
rule produces values in that range. From the theoretical 
point of view it seems easy to obtain similar results using 
a fixed connectivity and changing the bias, or modifying 
both of them simultaneously. 

Much future work is open. We think that the order- 
control coupling parameter proposed in [Tl| ] also works 
for RBN. And vice-versa: the proposed method in this 
article works in this type of binary neural nets. In or- 
der to consider biological implications the automata (to 
be considered as genes) dynamics and the disconnection 
rule (to be thought of as mutations) should have differ- 
ent characteristic times. We should also include a rule 
that allows for a possible connection growth, and not 
simply consider connection loss as we have done in this 
article. We think that the point of view considered here 



can have further applications along these lines. In any 
case, we can already extrapolate our results to a different 
type of nets that show a clear transition between the or- 
dered and disordered phases with well defined order and 
control parameters. 



Acknowledgments: The authors would like to thank 
C. Molina-Paris, U. Bastolla, and S. C. Manrubia and 
the referee for their critical reading of the manuscript. 
This work has been supported by Centro de Astrobiologfa 
(CAB). 



[1] S.A. Kauffman, The Origins of Order: Self-Organization 
and Selection in Evolution, Oxford University Press, 
1993. 

[2] S.A. Kauffman, J. Theor. Biol. 22 (1969) 437. 
[3] A. Bhattacharjya and S. Liang, Phys. Rev. Lett. 77 
(1996) 1664. 

[4] U. Bastolla and G. Parisi, Physica D 115 (1996) 203-218; 
ibid. 219-233. 

[5] R. Albert and A. -L. Barabasi, Phys. Rev. Lett. 84 (1998) 
5660. 

[6] S.A. Kauffman, Scientific American 265 (1991) 78. 
[7] A. -L. Barabasi and R. Albert, Science, 286 (1999) 509. 
[8] S. N. Dorogovtsev and J. F. F. Mendes, Phys. Rev. E 62 
(2000) 1842. 

[9] S. Bornholdt and K. Sneppen, Phys. Rev. Lett. 81 (1998) 
236. 

[10] M. Paczuski, K.B. Bassler, and A. Corral, Phys. Rev. 

Lett. 84 (2000) 3185. 
[11] S. Bornholdt and T. Rohlf, Phys. Rev. Lett. 84 (2000) 

6114. 

[12] K.E. Kiirten, J. Phys. A 21 (1988) L615. 
[13] B. Luque and R. V. Sole, Phys. Rev. E 55 (1997) 257. 
[14] S. Maslov and Y.-C. Zhang Phys. Rev. Lett. 75 (1995) 
1550. 

[15] M. Paczuski, S. Maslov, and P. Bak, Phys. Rev. E 53 
(1996) 414. 

[16] D. Sornette and I. Dornic, Phys. Rev. E 54 (1996) 3334. 
[17] P. Grassberger and Y.-C. Zhang, Physica A 224 (1996) 
169. 

[18] D. Sornette, A. Johansen and I. Dornic, J. Phys. I. 
(France) 5 (1995) 325. 

[19] B. Derrida, in Fundamental Problems in Statistical Me- 
chanics VII, edited by H. van Beijeren (North-Holland, 
Amsterdam, 1990). 

[20] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. 
59 (1987) 381. 

[21] For reviews see: P. Bak, How Nature Work, Copernicus, 
New York 1996. H. J. Jensen, Self-organized Criticalitity, 
Cambridge University Press, 1998. 

[22] B. Luque and R. V. Sole, Physica A 284 (2000) 33. 



8 



